{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demo Foreground subtractor\n",
    "\n",
    "Example notebook for subtraction of lens light to reveal faint lensed source features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from StandAloneRingFinder import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from astropy.io import fits as pyfits\n",
    "import pylab as plt\n",
    "%matplotlib inline\n",
    "import colorImage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(13.0, 13.0)\n",
      "(16.0, 13.0)\n",
      "(16.0, 13.0)\n",
      "(16.0, 13.0)\n"
     ]
    }
   ],
   "source": [
    "pref=1\n",
    "ddir=\"example\"#Change to your data directory\n",
    "imgg=pyfits.open(\"%s/clip_%s_g.fits\"%(ddir,pref))[0].data.copy()#[a:-a,a:-a]\n",
    "imgr=pyfits.open(\"%s/clip_%s_r.fits\"%(ddir,pref))[0].data.copy()#[a:-a,a:-a]\n",
    "imgi=pyfits.open(\"%s/clip_%s_i.fits\"%(ddir,pref))[0].data.copy()#[a:-a,a:-a]\n",
    "imgz=pyfits.open(\"%s/clip_%s_z.fits\"%(ddir,pref))[0].data.copy()#[a:-a,a:-a]\n",
    "\n",
    "sigg=pyfits.open(\"%s/clip_%s_weight_g.fits\"%(ddir,pref))[0].data.copy()#[a:-a,a:-a]**-0.5\n",
    "sigr=pyfits.open(\"%s/clip_%s_weight_r.fits\"%(ddir,pref))[0].data.copy()#[a:-a,a:-a]**-0.5\n",
    "sigi=pyfits.open(\"%s/clip_%s_weight_i.fits\"%(ddir,pref))[0].data.copy()#[a:-a,a:-a]**-0.5\n",
    "sigz=pyfits.open(\"%s/clip_%s_weight_z.fits\"%(ddir,pref))[0].data.copy()#[a:-a,a:-a]**-0.5\n",
    "\n",
    "psfg=pyfits.open(\"%s/clip_%s_psf_g.fits\"%(ddir,pref))[0].data.copy()\n",
    "psfr=pyfits.open(\"%s/clip_%s_psf_r.fits\"%(ddir,pref))[0].data.copy()\n",
    "psfi=pyfits.open(\"%s/clip_%s_psf_i.fits\"%(ddir,pref))[0].data.copy()\n",
    "psfz=pyfits.open(\"%s/clip_%s_psf_z.fits\"%(ddir,pref))[0].data.copy()\n",
    "\n",
    "\n",
    "color = colorImage.ColorImage()\n",
    "colorimage = color.createModel(imgg,imgr,imgi)\n",
    "\n",
    "psfmode=\"dont\"\n",
    "RF=RingFinder(imgg,imgi,sigg,sigi,psfg,psfi,0.265,1e12,1e12,visualize=False,psfmode=psfmode)\n",
    "RFgz=RingFinder(imgg,imgz,sigg,sigi,psfg,psfz,0.265,1e12,1e12,visualize=False,psfmode=psfmode)\n",
    "RFrz=RingFinder(imgr,imgz,sigr,sigz,psfr,psfz,0.265,1e12,1e12,visualize=False,psfmode=psfmode)\n",
    "RFiz=RingFinder(imgi,imgz,sigi,sigz,psfi,psfz,0.265,1e12,1e12,visualize=False,psfmode=psfmode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2IAAAEQCAYAAADS2sdxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt41eWd7/3vIidyJCSQlUACiUCQQ0yiFESLRhEUO0Uc\nZqxOx81W2un0mc60W/dMnT/2M+mefXWw+6kz2rpnujvWTe2M1T2OSFtB8BBFRBBJOAUIhwRyPhLI\n+bieP7xwPGDuD8nKLyvh/bquXFeRd+/fL2v91r3WTbLW7QsEAgYAAAAA8M6ksT4BAAAAALjasBAD\nAAAAAI+xEAMAAAAAj7EQAwAAAACPsRADAAAAAI+Fj8agPp+Pj2IEJqhAIOAb63MYCeYnYGJibgIQ\nioaam0ZlIWZm9sgjj9gTTzzh7HzzUqXxvlmRKHU/7zsudZlTrpW65gt1n/pzj3VblE3+XBcfF+Ec\nq7q9UTrmWPnSIq07e+oWZzM9sV4aqyGmQeoay89L3X/OC5O6f6+Ik7qLrRekThKR9ek/D5w3C5v6\n+S62RRouMlI7t17tJsblxEdJ2W1tsc7mbdPu10GpwuVp87pZtVRdf+N0qSvee+bj/x0ImPku85Rb\nsmm1NFbu93dIXfClfObPHWZ2mes6us091GCXdsgeLbualAU2Wpu99/Gff1bYaN8q/Px1GCXeeO/Z\nTVLXY9pct8z2CmNFfurPvyistYcL0z7XtVu8dMxo066n5+xBqcuzEqlbaKWf+vPmwkrbUJjxue41\nu9M51n+1/086ZqktlLoui5a6ts/cxv+3sMz+sDD7c12x5Uvj3WGvO5sMq5TGqrBMqVtV8+6n/lz4\nY7PCRz/fnZ2hzdc7hPvLzGy1vSZ1rut4se/MkH8/7F9N3L59+13XXnvt8Xnz5p18/PHHvz/ccQAg\nmJibAIQq5icAnzSshdjAwEDYd77znZ9u3779rtLS0oXPP//8A8eOHVsQ7JMDgCvB3AQgVDE/Afis\nYS3E9u3bt3Tu3LmnMjMzKyIiIvruv//+X7/yyiv3BPvkQlHY6P02J642vs//iitG5mqem4Dgcv+6\nPa7McOanGwpivDq9UZNfoL0VINTlFiSM9SmM2MKC5LE+hRErWD7WZxBcw1pVVFdXz8zIyPj4l0DT\n09Or9u7du+yTzZ49e0Z6biEpnIUYgmWS9jve0ClzEzCaLvf+sPEp0p3giijz008KD1iPffR+8hsK\nYmxJgfv9p6Euv0B7L1ioyyuYMtanMGKLJsJCTHvb45j5oKjLPijqlvthrSqUT/ZZvnz5hF2MAQhN\nfOoYgFClzE9/Xni9tZn+Ig5AaPlSQbR9qeA//qH9f/1g6A+bG9avJs6cObO6srLy44+OqayszEhP\nT68azlgAECzMTQBCFfMTgM8a1kJsyZIl+0+ePDmvoqIis7e3N/KFF1742tq1a7cG++QA4EowNwEI\nVcxPAD5rWL+aGB4e3v/Tn/70O3feeedrAwMDYRs3bnxmwYIFx4J9cgBwJZibAIQq5icAnzXsT55Y\ns2bNtjVr1mwb6QnMPtMpdUcHtM1Q1TeILJry2c0rL+/9C01St7Dd3fnNL41VbNpmyH87c67UPVd9\nSuqOndIuh4U97s0V99WfkMa65qvrpK6xfIvUne8akLqgbtRsZjYz3ZlE1mjXUm+rsGGqmfVK1dUn\nWHOTmZm1aZumnhY2V1U3ap7v0zYlPhHQNq9XqU8GyiNs0+z50lg/PavNE9r2oGbTxE1um0x7jJWF\na5/zEggMvWGnmdnzb+6Wxgq6ZG2T0wzxUxIrm8XNmnFZrvnpnM22Wmt1jqNu6Nxk06Su4XMbel9e\njh2WOsU7dot4zENSlyjcbmZmrZYodcpGzWZmS2y/szltc6Sx3rDbpe679qTU/da+KnVtpn2ypbLh\ndLw4v2ZahdQdmKHt8JBp5VJ3Xrz/ByxM6jptZJ9sOuwNnQEAAAAAw8NCDAAAAAA8xkIMAAAAADzG\nQgwAAAAAPMZCDAAAAAA8xkIMAAAAADzGQgwAAAAAPMZCDAAAAAA8xkIMAAAAADwWPloDp6Yk26L5\n1zi7oyfOSONViMeNi9W6V6vekbpB8bjlWe5d6ZvbLkpjhXdMlboDMyukLj05WuoaAj1St+/wB87G\nr21Ibi2NR7VQVNIb1OHsoTlat7e8ydlUzOyWxppWNUPqaqxG6iZHTpa67l7t/CaCuwu+JHWvFrmv\ndTOzc0IzSfxnrxODx7VQdNvya6XurT3BO+72mBNSVxm0I36kycqDOl77u68Hbaz32qLEskOqvrNI\nG62pf4rUvR9Zrw3YrGUYnkQ7b4Mm3hdBFG1dUrfN1jib9faSNNYy2yt1qj+256Tut/ZVqZshPsf2\nm/sFT7llSmNFWp/UVVm61M2xU1LXadrrxArLcjYHLVca63r7UOoSrVXq2ixe6jLEZ57TNlfqmizZ\nURwZ8m/5iRgAAAAAeIyFGAAAAAB4jIUYAAAAAHiMhRgAAAAAeIyFGAAAAAB4jIUYAAAAAHiMhRgA\nAAAAeIyFGAAAAAB4jIUYAAAAAHgsfLQGnjngs+v73Ou8o0E+7jX5c6Tu0Lung3rcc6mznE1v+X5x\nNG2H+yb/DKl75zfa7vCx102Vurm/v8bZdJVrO6bXv39C6v4o+6tSt3dBi9TNadwtdamDk6WudLDb\nHVVJQ1mnafdXwuQVUnexe5fUzU5yN2e1mzfkzb/YJ3WvBvGY3/iTO6Tuf//T60E8qtlb1dp8EtRj\nHgvygHFiN6dA6w4WDfNELu8vb/5jZ/M/j26RxkoQjxnW7JO6X9edEkcMpmSxax7VsxiPGm261ViG\nszth86XxbrL3pC7KeqTukOVInWKZ7ZW6GtNe61ywRKmbY9pj4k1bKXV/bCXOZkmP9vqvNGqh1NWK\nt0mitUpdgb0tdVvM/Tqx1yKlsdT7NUZ8TRxm/VL3lQHtmf1EmPYYe8G+5iheGfJv+YkYAAAAAHiM\nhRgAAAAAeIyFGAAAAAB4jIUYAAAAAHiMhRgAAAAAeIyFGAAAAAB4jIUYAAAAAHiMhRgAAAAAeIyF\nGAAAAAB4LHy0Bk46H2FZ56KDNt7kxFip8/fVSN3KWdru5Q3n3LuIm5mlVVY4m93h+dJYa7MPSd3s\nuMlSV5k3Reoeu3Cj1B0dcN+vNWkZ0lizK/KkblpWndRNSdYu6T3aRu32r1F+qStYOOBsqrqapLHq\n47X74U/LtdvkR91SZudatG4iqDhwLqjjxSbFOZtZvkZprD+5uUDq/vfuPVKXea5a6gaSvix17S3v\nOps2aSSz/iity2r3SV35wbPikTU3mDZn+xOF552L7dJYMVJl9qT28NeFid2Ae078ujjUv4jd1aTL\noq3N4p1dsmnPJ/Hio/Gw5UjdQit1NjHWKY21w1ZL3Vfst1L3D/Y9qVNvkzW2TerKLNvZTIvS7q9O\ncQaINu1FTKlpr3WTrVnqlGtznW2RxqqxNKnrNe2JYoEdk7pnwx6SukRrlbocOyx1X4SfiAEAAACA\nx1iIAQAAAIDHWIgBAAAAgMdYiAEAAACAx1iIAQAAAIDHWIgBAAAAgMdYiAEAAACAx1iIAQAAAIDH\nWIgBAAAAgMfCR2vgnYO19kT/yHab/qRbJkVK3aHT56WuM1AjdYujF0hdatigs0nIG5DG+lLkDKlL\naOmQulm3/57URQUSpe7anjPuseYlSWPFzpcyizrbL3WDAe24h9dPk7qBmmVSd3P6DmeTcG6xNNak\n9w9J3Y8GW6ROFQjqaKHtZQvubdfb0u5snvjXg9JYA5FTxaPOkaoYuyB1xxObpG6WcNPNk0Yy2zd9\nptSVV9WLI5Zr2SyflH1YrT2HHfidMj9px6yLER+JnVq36AZtuMqjCVLXM+C+L/5FOyQuI9mard/c\nr0+irVMar9QWSl2bxUtdjHU5mzhrk8aaIs5NLZYsdb0WJXVnLEXq/pM9J3XKbfye3SSNtcBKpa7L\nYqSuU+zm2GmpW2tbnc2ycu01zP/M+jOpWyjeJm/aSqmLtB6pazXtNbE63hcZ9kIsMzOzIiEh4WJY\nWNhARERE3759+5aO6EwAIEiYnwCEIuYmAJ807IWYz+cLFBUVFSQlJQX3n5YBYISYnwCEIuYmAJ80\noveIBQIB7XctAMBjzE8AQhFzE4BLRvQTsTvuuOP1sLCwgW9961s/++Y3v/nzT/79nj17Rn52ADAM\nrvkJAMaCa256tvAd67A6MzNbWJBsiwq090cBCA1lRXV2sqhO7oe9ENu9e/fNaWlptY2NjdNXrVq1\n89prrz2+YsWKXZf+fvny5SzGAIwJ1/wEAGPBNTc9VHiL1VvJWJ4igBHILki17ILUj//8ux8M/UFd\nw/7VxLS0tFozs+nTpzfee++9L/OGUwChgvkJQChibgLwScNaiHV2dsa0tbXFm5l1dHTE7tixY3VO\nTk7wPqseAIaJ+QlAKGJuAvBZw/rVxPr6ev+99977splZf39/+Ne//vV/Wb16tXsjJQAYZcxPAEIR\ncxOAzxrWQiwrK6u8pKQkL9gnAwAjxfwEIBQxNwH4rGF/WEewRIu7fu9oOS91EZO14/Z1t0rdofCX\npK5h9mxn092u7az+cs5tUndzRIXUfVfcRb7u+nSps7b5zmS2f0Aaqn+a9tuxk1O17paqNqnL6p0u\ndTsXR0pder/7/PrPHpHGOjNVu7+sWduG5p54bbhXtJvuqhJuUVLXZz3O5kKXdsyBC9pcF2VaN1hw\njdYVHZe6iqzb3U3jAWmslKo+qRv48s1S13ziotRZRKfWDWi3cSDOPf/727U5sbtT+7StC9YrdRG9\nTVI30K3ddu4rXbdGfL7e1h3Eg4a4JGuxSVYTtPGSTbv/Vb+xrzqbKPEqWWlvSN28hiqp29TxA6l7\nJ+tLUnfYcqRuVcO7zub2niB/eJ34mLgj1n1uZmbnZmiviSosy9l0ztB2Z2i1qVL3thVI3a1WJHXx\n4vP6q3aL1C2zvVL3RUa0jxgAAAAA4MqxEAMAAAAAj7EQAwAAAACPsRADAAAAAI+xEAMAAAAAj7EQ\nAwAAAACPsRADAAAAAI+xEAMAAAAAj7EQAwAAAACPhY/WwNfOvtbW3LjG2YWJu35HVmZL3ZTuMql7\n9pyUWXu/1k1pH3Q2N93UKI3l671W6pIy0qQuISNF6rrS0qUuYnKcs/G3Jktj+WIDUnd+dpPUZR6u\nlroZLdr3OnOqdj1NLv6ys3njnq3SWNUvtUidL13blf5UT4zUzexvczbVXRelsULfdKlatlS7L1Ja\n5zmbirKT0ljFUmXWI3b+/jCpK5+vXU89tQPuKOCTxppb4L7dzMzeUy+7GUla1zxHyhbcqs2xrTHu\n+enLDZ3SWOWVvVLXP6dO6lLKpkldiWlzrCI+O1XqtpVp38PVJMp6LMa6nN0uWyGNl2zNUrffbpC6\nAeFl4532mjRWp2nPTW+lLJe622r2SN008VoPM2GuMzN7Q2iOa0NZudiJr959eeJ4fyF2gn+L+gOp\nm2E1UnfCtNf/qoVWKnWJdl7qei1qJKfDT8QAAAAAwGssxAAAAADAYyzEAAAAAMBjLMQAAAAAwGMs\nxAAAAADAYyzEAAAAAMBjLMQAAAAAwGMsxAAAAADAY6O2oXNmR5fd1OzeDO2pQ+5NZM3MllyndeL+\ncLZyhtZNSlgsdfcluDcSjW3qlsa6+HXtm5jXc53U9V0zU+oW98RKXVxcpLOZnKJthNgToV2CqR1+\nqetcoG1K29mkba6adipR6urj3N9Ha+oqaaw5N2hb9e5/9R2pOypuSng1uca0zdVL9mvjTZ3q3qz5\nD5fGS2NNOafNdUXiXrjvvKttJP3ni7VNzp9KetsdNWdIY51M0DYbtgptY+3pKZOlbvJcbdPkxNla\nd+Gse2P6munabbL/Q2V3WDNraNc6UfI92v3f/EqVs2kL8kbNPmF/8ID7LhgXTto8qxa2a28zbT5p\nNe05bJq48fOAuTeIVzfCLbP5Uqdu/BwQX9ctGdAm9oR/7tMG7Bca9832kRSxOyN2RVo2u0Z7TozZ\n5N5svMIypbGSxY21/8AOS5163amUjdXNzLrE6/OL8BMxAAAAAPAYCzEAAAAA8BgLMQAAAADwGAsx\nAAAAAPAYCzEAAAAA8BgLMQAAAADwGAsxAAAAAPAYCzEAAAAA8BgLMQAAAADwWPhoDRyIqbHBpIPO\nbq64cXnxualSl2ADUjc163apm5T+htRNr3I3XctzpbFiUnKkLq56kdRNC++Wutj0XqlLupDkbCIm\n+6SxWuITpO5CUpTUTZrcKXVT4jOkrravTeriOw84m4yyVGmsi+m1UqfK1m5iK7sY1MOGtDuu0br+\n2X6p+/fDHc6mNna2NFbuPaulruhnT0iddtWZPVUqTGJmZlNjnUls3A3SUCldAanzze2XusRW7Smt\nq1WbE8+I80lW/gJns3vvKWms1DuWS13doRKps4ZGKUs/qV0pzeHCdaLdXbKAdplMCOWWZcfN/fzp\ntwZpvBtsv9S1WbzULbc9zuZNWymNVW6ZUrfEPpS624+7z83MzH6jZZYmdu6XRGYnxbFEH76sdb8V\nx/sj7SWbzVvZ7myWr9Luhy12j9TNtdNS12UxUvdv9gdSN8NqpG6/LXEUQ98L/EQMAAAAADzGQgwA\nAAAAPMZCDAAAAAA8xkIMAAAAADzGQgwAAAAAPMZCDAAAAAA8xkIMAAAAADzGQgwAAAAAPMZCDAAA\nAAA8Fj5aA5f0TLcX2rKd3WDFcWm8efMXSF1fWIfU3XEuWuqO5Wo30ZlJmc4mMzBdGmtBjXZuMYlS\nZkkJ2m7jiZPatAFThNs4TPsekrq029eXoP2bQVR8t9Ql1CdIXbt43PZ5F5zN1LpKaaxAQ6rUfTNz\nstRt69MulKzkOmdTXi4NFfJ+2TxT6qadqZa6NnPfZ795t1Ya66Hi01JnYVrWMKB1k6K0eWKw2d11\nJLZIY7W2R0ldb32V1JW1aV3KtBSpa/iwR+rip7ofY9ERF6WxUrq1+6FuUq/U2Qwta63br4X97iRX\nG8kOit3V5DZ7y663953dPlsqjTfNmqUuxrqk7ozNcTZtFi+Ntcz2SZ3f6qXu4twIqUuI7JM6a9Cy\n6q8nOZuZydqcaGValu8+pJmZ/VY87Lz5Wmft7qTPtPuh1aZK3Xt2k9QV2FtSF2na3Jljh6VOvd6/\nyJCvMh9++OFf+P3++pycnI/PpqWlJWnVqlU7s7Ozy1avXr2jtbVVXA4AQHAwNwEIVcxPAFRDLsQe\neuihZ7dv337XJ//bpk2bHlu1atXOsrKy7JUrV76xadOmx0b3FAHg05ibAIQq5icAqiEXYitWrNg1\nderU85/8b1u3bl27YcOGzWZmGzZs2Lxly5Z1o3mCAPBZzE0AQhXzEwDVFb9HrL6+3u/3++vNzPx+\nf319fb3/ct3OnTussfGj95zExMRZbGzciE4UgPe6usy6tbfdjTl1bjIz6+/+j/frTAqPsknh2nuV\nAGA41Pnp6cL91mMfvXbKL4iz6wt47QSMJxVF56yiSPtMALMRfliHz+cL+Hy+wOX+btWq1VZTo73B\nEkBoio7+6OuS1taxO5crMdTcZGYWPln7sBYACLah5qc/K1xiF5RPRAEQkjILZllmwayP/1z0g91D\n9lf88fV+v7++rq4u1cystrY2LSUlRfxcGQAYPcxNAEIV8xOAy7nihdjatWu3bt68eYOZ2ebNmzes\nW7duS/BPCwCuDHMTgFDF/ATgcoZciD3wwAPP33TTTe+dOHFifkZGRuWzzz770GOPPbZp586dq7Kz\ns8vefPPN2x977LFNXp0sAJgxNwEIXcxPAFRDvkfs+eeff+By//3111+/Y3ROBwDcmJsAhCrmJwCq\nEX1Yx1AuRPVbVZz749YuWo80Xnz0Malb2DfLHZlZ8srnpW7JrilSl3h7lrPxx6dJY01J0m6TpNgv\n/CyCT0mIjHZHZma9nVo3pdbdDGq3mw1ob0qe6tM+XCGirU/qTvu6pM53Qfw+jsU4k7rYImmogbIM\nqUtdoF3rVdnip249Wad1E0D3hWqpqxLHm+5333aN4mcXlSS8InVp4genpPgXSd2i62dL3WvHDjub\njnbt8RUfiJW65jPq22m0ubOhUftEq8Qp2mPsVKl7Hksc0ObhQ401Umf1wf1Ah7NBHGvazDAtDB/Q\numCeXIiLsU4L2EWpU5yw+VIXKT52XrT7nE2meIflWYnUxVub1JWHZ0pd7qyTUmdHtOxV393O5mfL\n/kQa685lO6Tu+/c/LnV/c732mshu1DITPmA4zbQ5LNHOuyMzO2VzpC7ZmoN6XPX6bLN4qfsiV/we\nMQAAAADAyLAQAwAAAACPsRADAAAAAI+xEAMAAAAAj7EQAwAAAACPsRADAAAAAI+xEAMAAAAAj7EQ\nAwAAAACPsRADAAAAAI+Fj9bAadM6LWeue/fqWeHR0niviDthz5ycKXVnLkiZ9S9OkbrkiF5n4wtv\nkMaaHOuTutjELqkzE3dWj9RuY+sRbrzYP9TG8v1O6zp6pCyu1X0/mJn19Wnfa1N4rdRNb3XfFzFn\nV0tjfa1yr9T9Q0a21GW99K7URSW4m+MXpaFC3wwtS2/RuqoOd7M4O10aa09ZldRNmSZldrD+qBae\n1R4T12RGOZtztdrjtbbskNR1mzaez6fNnf39/VIXFhamdfXu47bWifNwdLfWJQ1IWaR4DV+XnCF1\n+5srnc0b1doxVVOFRnz2Cnmnba7VmntCqTe/NF6Gue8vM7N4a5O6TDvrbLLthDTWgGmPrxOmPddp\nj36z3LyTUlexTnv99639v3Q2gaXanHPAt0zqTg/Mkbpfv/SQ1FmclimrhkrT5pIG8RousLel7rf2\nFan7a9skda/a3VL3mt3pKLYN+bf8RAwAAAAAPMZCDAAAAAA8xkIMAAAAADzGQgwAAAAAPMZCDAAA\nAAA8xkIMAAAAADzGQgwAAAAAPMZCDAAAAAA8xkIMAAAAADwm7JE9PL1d7dbZ2uDsJkfMl8aLqtb2\nTC+ZP0vqvrozQuoCa49rXdIHzqbh9O9JY81NGJC67pgWqYvtFrdMj+vUuqg+d9P8hDZW+xSta2vW\nupYeKZvSdFHqZn6gXSdHUgLOZvK/1UtjnZ+eLnXdAW288irtNrl7kbs5flQaKuTl5Wpdz75ILWzu\ndSZHLFUaaoVVSd2uJimz6Rlad7BUfIyVCtfTtCRtrBZtDlMFAu7HoZlZeHhwn/oi6s46m9kxydJY\n1d39Ure4b47UJa7Qrqd3d1VKnXIV18WI/8bbqT3XDWqjTQjNlmQ1lubsZliNNN4S2y91qVYrdYft\nOmfznt0sjVVl2uS03l6SuhLLk7qTWdpzbJEVSF3gQ+H1aUCbcwL236Xuxf/1/0pd5/8TI3XrbIvU\nZZh7nggz7XGtjGVmVis8HszMsu2k1LWY9vz0qt0tdadMm4u/CD8RAwAAAACPsRADAAAAAI+xEAMA\nAAAAj7EQAwAAAACPsRADAAAAAI+xEAMAAAAAj7EQAwAAAACPsRADAAAAAI+N3obOleHW8UGUs2u/\npk4aL8eETYTNbMZ+bYPA0/naZojhJbFSlzHzgrNJG3Rv+mlmdr5S2+R6IFzbgDk2WdtczyLc34OZ\nmbWEuZuAtimpXRQ3dO3R7oemKm3rz6RabZPTwCTt/N5rOe9siuO0DZgrL2obMKbGnZG6/Fv8Urdj\nr3Z+E0HJNq27Kcu9UbOZmQl7ISeWaXPOLu2I8r+itWp7ZtpNk7Sng/cGhQ2dm7SNlaf44qXuQqBN\n6lTqxs+qbqFp7NSew9bEa5vId7Zpm5d+eES7hrVt382UZ+yESO17vag9hVmPsF+uBfcuHTNx1m5T\nrdXZtVqiNN5bdpvU5dhhqRsw9/N/ppVLY6nfw+u2UuoWWqnUtZk270wx8TXRncL1/pz2uLaDf6N1\nc7QLfqntk7oU057/X7M7pU6hbkoeY9pEoV53L9jXpG6F+Gx8g2PT9AOO/z8/EQMAAAAAj7EQAwAA\nAACPsRADAAAAAI+xEAMAAAAAj7EQAwAAAACPsRADAAAAAI+xEAMAAAAAj7EQAwAAAACPsRADAAAA\nAI+Fj9bAUe3JFtd4jbO7kNEsjRdI0HZgj0k+LXU1dV1SN2kgQ+paMqOcTXLkB9JY/We1Xd+tI1LK\nBpc3SN2MFm0864t1N9FN2liTtN3h26rDpK71YqvUlZzXjtsXp+3UHn2qz9l0pKZKY8XUn5G6g3HT\npe5Uk/txaGY2p++is/nQtMfNRPGedvdL0uPqpK61XRvvrkTtfn01T7ueztV1SN21J6Y4mwrt4WUX\nfG1aKI6XoGXmnq0/0ih2ii5zP77MzLaIN0nqVK1rPq91wXRRu5TMTHteXxLtvgDe7bygHjSkzbBa\nm2zu1zEVlimNl2jac2KpLZC6SOtxNnPslDTWgPgS9HVbKXXTTHs9GSM+jx0Tb5PbM19zNoHZPmms\ntxZ8Rer8N2hPTvWWInW77BapW2Z7nc0hy5HG6rQYqVtopVJXadrr9XW2RerSrVLqttndUvdFhvyJ\n2MMPP/wLv99fn5OTc/jSfyssLCxMT0+vys/PL87Pzy/evn37XSM6AwAYBuYnAKGIuQmAasiF2EMP\nPfTsZycLn88XeOSRR54oLi7OLy4uzr/rrru2j+4pAsDnMT8BCEXMTQBUQy7EVqxYsWvq1Kmf+8WG\nQCCg/YwVAEYJ8xOAUMTcBEA1rA/r+MlPfvLnubm5Bzdu3PhMa2vrZX/J+7U9O63WGq3WGq3N5F8Y\nBxBC2gIDVhPo+/hrPFDmJwDj24WBfjvX1/3x13igzE3/VPi+bS6stM2FlVZSNDHe9wZcTU4W1dqr\nhQc+/nK54oXYt7/97X8sLy/PKikpyUtLS6t99NFHf3y57s7lqyzNpluaTbd4Ez7cAUDIifeF2Qxf\nxMdfoU6dnwCMb1PCwm1WxOSPv0KdOjf9aeGNtqEwwzYUZlhegfuDcQCElnkFaXZ34fUff7lc8UIs\nJSWlwefzBXw+X+Ab3/jGP+/bt2/psM4UAIKM+QlAKGJuAnA5V7wQq62tTbv0v19++eV7P/mpQAAw\nlpifAIRHmJjnAAAUo0lEQVQi5iYAlzPkJg4PPPDA82+//fatTU1N0zIyMip/8IMf/E1RUVFBSUlJ\nns/nC2RlZZX/7Gc/+5ZXJwsAlzA/AQhFzE0AVEMuxJ5//vkHPvvfHn744V+M3ukAgIb5CUAoYm4C\noNK2NR+GKqu3A3bM2eWeiZTGO9jaIHWBpFlS9/vx2g7snT5tV/q6F/zOpmNpmDTW8WXvSd3sw6el\nzn8qzR2Z2YyqBKnrvqbG2Uzu0nau72jTLsH4jhNSt79Nu/8rB7U3QZdXzJW6+He3Opv5/dq19MGU\nQam7vTJe6rac2iN1XTZdqq4m7kf1R+qF5ki7dr9miMd8q/WMFhZpWYT28LfjucInuZVoY5l2k9h8\ncTj183nV21j93X1l1lE//64qRnkcmtV1R0tdtJ2Tululykza/Er+gFVtTny3Ux1v/Guw6VYrXKGd\nFiONt9Jel7p9tkzqdtktzqZMfMRmm/a8Hm4DUndCPO4p057Xw6xf6jKs0tk0+LRnk6S8aqlbYh9K\n3a96HpS6v4r6kdSlmPu1+Hp7SRrrRfua1P3O7pa6GVYrdW2mvXZS7/9ltlfqvsiwPr4eAAAAADB8\nLMQAAAAAwGMsxAAAAADAYyzEAAAAAMBjLMQAAAAAwGMsxAAAAADAYyzEAAAAAMBjLMQAAAAAwGMs\nxAAAAADAY+FjfQKL+mZJ3cGUKVK3u0HbbbxZ2+Tc+i8kSl3brk5nszQyShqrY+YBqTtard12PQFt\nd/AzDdp4i7vanU3ngHbMurMDUtebJWVWddq967uZ2ey2rVL31rRbpa5hxqCzWVrnbszMstOkzP5t\nT5vUrZimjberrdkd9WhjTRTa1Wmm/ZuWdv9XysdUaf/edl2Hdn7h0e7mpDSSrsfipO6iuecmM7Ma\n8bj1YhfMsdI7G6XO/YzzkW6xu3Wx9hxb33jB2RSL32x+jNYVDwjXcI92/Ya6Lou2Not3djHiFeDv\n0K6nnljt9ckcO+Vs8q1EGivauqSuzOZLXYnlSV2OHZa6SKkyW2lvOpt37cvSWHPC3LevmVmWVUhd\nRpT2jKJeTz+yv3I2/83+VhprrnAtmZkVi/drsyVLXZPYtYvPO/fZi1L3RfiJGAAAAAB4jIUYAAAA\nAHiMhRgAAAAAeIyFGAAAAAB4jIUYAAAAAHiMhRgAAAAAeIyFGAAAAAB4jIUYAAAAAHiMhRgAAAAA\neCx8rE/g9Wu0ndXtsLYDd0A8bmf7TKlb2hEmdd2BWGfzTKI0lC1rrpe6HF+E1NWcvSh15WEdUlfX\n5b4v0iqipLFOZw9IXcPk81I3pVfrDh3R7tfTt74ldfnn053NL6O12+S2WU1SN7NRu9oP1QxKnfWI\n3VVEuyfMzLy/7eabdg23iNP81uQeqbtvobtZHdD+je/p97XbrcLapU6lzYhmPvEZst49/VvUBW2s\nWi2TzRf/vfWHbdr89Mg3M51N8f+okMYq65cys96rZ25aavssx952dpU2SxqvJDZX6g5bjtTdYa87\nmwrLlMZKtFapu8e2SN1v7atS12kxUpdszVI3y845mxlWI42lnlu/OP9nWrnUVVqG1C22I85GvZZU\n99mLUve2FUjdfDshdedNe9G+35Y4irND/i0/EQMAAAAAj7EQAwAAAACPsRADAAAAAI+xEAMAAAAA\nj7EQAwAAAACPsRADAAAAAI+xEAMAAAAAj7EQAwAAAACPsRADAAAAAI+Fj/UJNLQ0jMlxz1VWS93v\nLdTGezPO3SRkaLuN/58Dk6Xu/k6tm3Z0r9Rds2jo3b8vaWty7zZ+LKlHGmvmK71SV5J/UuqaW7U7\nrLMzV+qmxTwvdR+Yzz1Wc5c01o7SAam7sUPK7HSn1uEy4sWubVTP4rIqTbtOOsXOLmjZi6+5m/88\n4xptMDsldppZYlcvdn39WhcQbrt58dp8Hd52o9SVTC2SuuPn3XOTmVlmkzZRPLkteM/ZHdr0f1Wp\nt1RrsSxnl2xN0ninba7UrbeXpO6ILXY2u+0maazJpr1OuM9elLo7TZiczKzZkqXuHbtF6iotw9n4\nxVlHvV9rbYbUqXLssNQttX3O5jf2VWmsZaa9Nq2ydKlTb7tjtkDqNtj/kboSy5e6L8JPxAAAAADA\nYyzEAAAAAMBjLMQAAAAAwGMsxAAAAADAYyzEAAAAAMBjLMQAAAAAwGMsxAAAAADAYyzEAAAAAMBj\nLMQAAAAAwGPhY30CVt03JocNiN3TpVoXPyPM2Sx6tkIaq6G1W+qem67tIn6hUcrsm811UpcwaZqz\nOXzwiDSW72KU1K2IzpS6k4dOSt2/mHbHZu9Olbrw2Hpnk9jjvkbMzE63ag/Li1MHpc5qxQ6f1zbW\nJ/DFOoM83rxen9SdrHLPns9WnZLG0o5oFi9258ROPa76PLFGaPa1afN6sxVpBz2vZWYDUtXR0S51\nrR+6G3+yNJTVN2sdPu+0zZW6RGuVukHx3+UzrcLZZNsJaaxamyF1Jyxb6hbaMambbg1S95bdJnWn\nbY6zybFD0lg14m2SbNqDZ8C01x2ltlDqeizS2SwUX19VWobU9QrHNDPLscPieNrrzgFxiRQzwmfj\nIR95lZWVGbfddttbixYtOrp48eIjTz311F+YmbW0tCStWrVqZ3Z2dtnq1at3tLa2Jo7oLADgCjA3\nAQhVzE8AVEMuxCIiIvr+/u///r8cPXp00fvvv3/j008//WfHjh1bsGnTpsdWrVq1s6ysLHvlypVv\nbNq06TGvThgAmJsAhCrmJwCqIX/ulpqaWpeamlpnZhYXF9e+YMGCY9XV1TO3bt269u23377VzGzD\nhg2bCwoKij47oezZs2f0zhrAVW0kcxMAjKaRzE//VPi+dVm1mZnlFiRYXsEUz88fwPAdLWq20iL9\n967l94hVVFRkFhcX5y9btmxvfX293+/315uZ+f3++vr6ev9n++XLl7MYAzDqrnRuAgCvXOn89KeF\nN1qL7fP+RAEExaKCZFtU8B9vlP2/Pygbspfendne3h63fv36l5588snvxsfHf+rt6z6fL+Dz+dT3\nNANA0DA3AQhVzE8AXJwLsb6+voj169e/9OCDDz63bt26LWYf/UtOXV1dqplZbW1tWkpKivYRNAAQ\nJMxNAEIV8xMAxZALsUAg4Nu4ceMzCxcuLP3e9773D5f++9q1a7du3rx5g5nZ5s2bN1yaZADAC8xN\nAEIV8xMA1ZDvEdu9e/fNv/rVr/74uuuuO5Sfn19sZvZ3f/d3f/3YY49tuu+++1585plnNmZmZla8\n+OKL93lzugDA3AQgdDE/AVD5AoHg/4qyz+cLPPLof7Unnvh7dxzQNpsMtuXiRn177B2pi4pxf+5J\n3kJtl8uTUxOkLqvyuNTNEffMPtivdbNmuTeSjK/TNnT9d3HD3Jmr86Wu9lfFUjco7nGcdk261KVO\nde+uWvxhh3ZQUVbqEqmbn35R6rbvH/oNpZcEAgF1P9yQ5PP5AjZJ2+TSBr2fn9Yk3C11O9p3St3A\noDgBqFKEz3hq0CaTVe694c3M7ID4LWRmpknd5Kpaqdut7YVrliE8JCrG6O1AsWIX1OlptlSt0PaQ\ntV2lZ6VuIsxNjYHftx5hU+/X7Q5pzDBxQ+9W07Y0WyG8JpombjasbiKsnpu6kfSr9pWgHjfZmpxN\nlrARtpnZVlsrdV+zX0tdulVJ3Zu2UuqazD1p32mvSWMpG2Gb6dfwTfae1KkbRKvXp2vj73/0PTrk\n3KRtpQ4AAAAACBoWYgAAAADgMRZiAAAAAOAxFmIAAAAA4DEWYgAAAADgMRZiAAAAAOAxFmIAAAAA\n4DEWYgAAAADgMRZiAAAAAOCx8FEbOTBoFtB2w1ak2BSpa7ALUtdv+6Vu2rwkqetrC3M20b4+aawF\nDcel7obExVJX/P4RqdP2pDcbjLnobG7rnieNdaudlLq3f1ksdZMSo6TOWnukrPaMtit9T5z7/o+P\nnimNlRZZLXVlddo1XF4nZVeXweDNTcHWffEDqZuxeLrUVR6pGcnpfF5DvzPJEofa06R17eJ4PTUt\nUvfdmddLXd35A1J3uiLgbCZPi5HG6m7qlDqVr0PrYqLSpC4v1X2n7T57VhprV6mUXVX6LMJ6LNLZ\nxVubNF6MaddTppVL3S67xdksNO2OVb+HGpshdc/YRqnLMe01UbMlS92A8FI60nqlsW6y3VL3nP0n\nqZtp2uuJe+wVqfvQbnA26u0WZtrzcJ6VSN2/2h9JXaK1Sl2naXN2solPZF+An4gBAAAAgMdYiAEA\nAACAx1iIAQAAAIDHWIgBAAAAgMdYiAEAAACAx1iIAQAAAIDHWIgBAAAAgMdYiAEAAACAx1iIAQAA\nAIDH3NuBh4gGuxDU8T4Qd3S3k1q3QhlK/Baqu7SuNVnbHf6oNpzs5LEGZ3PK3I2Z2bovacd8IqB1\nP4ztkbombWN1WUu7e4f4yKhOaayy4F7qGGfeskYt1B7+do143Ivis0FTv7vpEsdqF8a6Eu2N2uP/\nh40HpG5BinbcJcJ0t79Je/wHmzh1Wkev9ly3+2zf8E8GTgct12qEfyPvtBhpvDwrlrpmmyZ1hy3H\n2VRahjRWm8VLnd/qpS7auqXunHh+d9prUldvfmej3iY32IdS957dLHUfmPYiK9q0F56rbYeziTJt\nHm61RKkbsDCp6xc79bETL64TVtg7Q/79/3D8//mJGAAAAAB4jIUYAAAAAHiMhRgAAAAAeIyFGAAA\nAAB4jIUYAAAAAHiMhRgAAAAAeIyFGAAAAAB4jIUYAAAAAHiMhRgAAAAAeCx8rE9grCzwad2xgNbN\nF3bq3tXVqQ0mOtoc1OGCSrzZbGey1vVkTJW6pn8/Lx5ZExmhdb3CQyk6sVsbq147JiYm8SFh6sP/\n3knaiD/uD96EUtcftKFGhTo/nc7VuuQ6ITosHlQ0Sfxn1IjJk6Wup6t9BGeDYJlnZZZm+51dqS2U\nxquwLKnrsmipy7QKZ7PQSqWx3rObpC7DKqVujp2Wug/thqAed8DCgtKYmVVZutTdbb+TulZLlLoG\n80tdpWU4G/V7nW8ngnZMM7N4a5M69X5daMek7g1b6SgODPm3/EQMAAAAADzGQgwAAAAAPMZCDAAA\nAAA8xkIMAAAAADzGQgzjUnNj31ifwoj19w6M9SkACLKeifKZF+onmmBU7S8K7od8jYUjRS1jfQpB\nsa9I+8CtUFZWpHyyUGibKNfTJSzEMC61NIX4x7IJ+nsHx/oUAARZb8dYnwEmkv1FXWN9CiN2pCi4\nn2Y8VibCQuzkhFiITYzr6RIWYgAAAADgMRZiAAAAAOAxXyAQ/F8E9/l8/HY5MEEFAgFxO/TQxPwE\nTEzMTQBC0VBz06gsxAAAAAAAX4xfTQQAAAAAj7EQAwAAAACPsRADAAAAAI95shDbvn37Xddee+3x\nefPmnXz88ce/78UxR0NmZmbFdddddyg/P7946dKl+8b6fBQPP/zwL/x+f31OTs7hS/+tpaUladWq\nVTuzs7PLVq9evaO1tTVxLM/R5XLfQ2FhYWF6enpVfn5+cX5+fvH27dvvGstzdKmsrMy47bbb3lq0\naNHRxYsXH3nqqaf+wmz83RcT0USYn5ibxg7zE0bLRJibzJifxgpz0zgRCARG9au/vz9szpw5p8rL\nyzN7e3sjcnNzS0pLSxeM9nFH4yszM7O8ubk5aazP40q+3nnnnRUHDhzIX7x48eFL/+0v//Ivf/T4\n44//VSAQsE2bNn3/+9///qaxPs8r/R4KCwv/5sc//vEjY31u6ldtbW1qcXFxXiAQsLa2trjs7OwT\npaWlC8bbfTHRvibK/MTcFFrfB/MTXyP9mihzUyDA/BRK3wNzU+h9jfpPxPbt27d07ty5pzIzMysi\nIiL67r///l+/8sor94z2cUdLYJx9PO6KFSt2TZ069VPbkG/dunXthg0bNpuZbdiwYfOWLVvWjc3Z\naS73PZiNr/siNTW1Li8vr8TMLC4urn3BggXHqqurZ463+2KimUjz03h6PJhNjLnJjPkJo2MizU1m\n4+vxYDYx5ifmpvFh1Bdi1dXVMzMyMiov/Tk9Pb2qurp65mgfdzT4fL7AHXfc8fqSJUv2//znP//m\nWJ/PcNXX1/v9fn+9mZnf76+vr6/3j/U5DcdPfvKTP8/NzT24cePGZ8bTj6UrKioyi4uL85ctW7Z3\notwX49VEmZ+Ym0IP8xNGYqLMTWbMT6GGuSm0jPpCbCJtULh79+6bi4uL87dt27bm6aef/rNdu3at\nGOtzGimfzxcYj/fRt7/97X8sLy/PKikpyUtLS6t99NFHfzzW56Rob2+PW79+/UtPPvnkd+Pj49s+\n+Xfj9b4YzybK7c3cFFqYnzBSE+m2Zn4KHcxNoWfUF2IzZ86srqyszLj058rKyoz09PSq0T7uaEhL\nS6s1M5s+fXrjvffe+/K+ffuWjvU5DYff76+vq6tLNTOrra1NS0lJaRjrc7pSKSkpDZcefN/4xjf+\neTzcF319fRHr169/6cEHH3xu3bp1W8wmxn0xnk2U+Ym5KbQwP2GkJsrcZMb8FEqYm0LPqC/ElixZ\nsv/kyZPzKioqMnt7eyNfeOGFr61du3braB832Do7O2Pa2trizcw6Ojpid+zYsfqTn0Qznqxdu3br\n5s2bN5iZbd68ecOlC3s8qa2tTbv0v19++eV7Q/2+CAQCvo0bNz6zcOHC0u9973v/cOm/T4T7Yjyb\nCPMTc1PoYX7CSE2EucmM+SnUMDeFIC8+EeTVV19dk52dfWLOnDmnfvjDH/71WH9CyXC+zpw5k5Wb\nm1uSm5tbsmjRoiPj5fu4//77n09LS6uJiIjoTU9Pr/zFL37xUHNzc9LKlStfnzdvXtmqVat2nD9/\nPnGsz/NKvodnnnnm4QcffPCXOTk5h6677rqD99xzz5a6ujr/WJ/nUF+7du36ss/nG8zNzS3Jy8sr\nzsvLK962bdtd4+2+mIhf431+Ym4Kre+D+YmvYH2N97kpEGB+CqXvgbkpNL98gcC4/bVKAAAAABiX\nPNnQGQAAAADwH1iIAQAAAIDHWIgBAAAAgMdYiAEAAACAx1iIAQAAAIDHWIgBAAAAgMf+f+n7wOYP\nK0r9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10eadec50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(15,9))\n",
    "ax=plt.subplot(1,1,1)\n",
    "size=colorimage.shape[0]\n",
    "\n",
    "#RF.Dshow[RF.fitmask]=numpy.nan\n",
    "plt.subplot(131)\n",
    "color.bMinusr = 0.8\n",
    "color.bMinusg = 0.4\n",
    "color.nonlin = 1.\n",
    "colored=color.createModel(imgg,imgr,imgi)\n",
    "plt.imshow(colored,interpolation=\"none\")\n",
    "ax.xaxis.set_visible(False)\n",
    "ax.yaxis.set_visible(False)\n",
    "plt.subplot(132)\n",
    "colorresid=color.colorize(RFgz.Dshow,RFrz.Dshow,RFiz.Dshow)\n",
    "d=0\n",
    "#plt.imshow(RF.Dshow,interpolation=\"none\")\n",
    "plt.imshow(colorresid,interpolation=\"none\")\n",
    "#plt.xlim([0,size-20])\n",
    "#plt.ylim([0,size-20])\n",
    "ax.xaxis.set_visible(False)\n",
    "ax.yaxis.set_visible(False)\n",
    "\n",
    "plt.subplot(133)\n",
    "plt.imshow(RF.Dshow,interpolation=\"none\")\n",
    "ax.xaxis.set_visible(False)\n",
    "ax.yaxis.set_visible(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
